// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtCracksDensityFilter.C
// ************************************************************************* //

#include <avtCracksDensityFilter.h>

#include <vtkCellData.h>
#include <vtkDataSet.h>
#include <vtkDataArray.h>
#include <vtkIntArray.h>

#include <DebugStream.h>
#include <avtCallback.h>

#include <string>

// ****************************************************************************
//  Method: avtCracksDensityFilter constructor
//
//  Programmer: kbonnell -- generated by xml2avt
//  Creation:   Mon Apr 30 11:40:50 PDT 2007
//
// ****************************************************************************

avtCracksDensityFilter::avtCracksDensityFilter()
{
}


// ****************************************************************************
//  Method: avtCracksDensityFilter destructor
//
//  Programmer: kbonnell -- generated by xml2avt
//  Creation:   Mon Apr 30 11:40:50 PDT 2007
//
//  Modifications:
//
// ****************************************************************************

avtCracksDensityFilter::~avtCracksDensityFilter()
{
}


// ****************************************************************************
//  Method:  avtCracksDensityFilter::Create
//
//  Programmer: kbonnell -- generated by xml2avt
//  Creation:   Mon Apr 30 11:40:50 PDT 2007
//
// ****************************************************************************

avtFilter *
avtCracksDensityFilter::Create()
{
    return new avtCracksDensityFilter();
}


// ****************************************************************************
//  Method:      avtCracksDensityFilter::SetAtts
//
//  Purpose:
//    Sets the state of the filter based on the attribute object.
//
//  Arguments:
//    a         The attributes to use.
//
//  Programmer: kbonnell -- generated by xml2avt
//  Creation:   Mon Apr 30 11:40:50 PDT 2007
//
// ****************************************************************************

void
avtCracksDensityFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const CracksClipperAttributes*)a;
}


// ****************************************************************************
//  Method: avtCracksDensityFilter::Equivalent
//
//  Purpose:
//    Returns true if creating a new avtCracksDensityFilter with the given
//    parameters would result in an equivalent avtCracksDensityFilter.
//
//  Programmer: kbonnell -- generated by xml2avt
//  Creation:   Mon Apr 30 11:40:50 PDT 2007
//
// ****************************************************************************

bool
avtCracksDensityFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(CracksClipperAttributes*)a);
}


// ****************************************************************************
//  Method: avtCracksDensityFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the CracksDensity filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: kbonnell -- generated by xml2avt
//  Creation:   Mon Apr 30 11:40:50 PDT 2007
//
//  Modifications:
//    Kathleen Bonnell, Wed Sep 29 09:03:04 PDT 2010
//    Changed output name to "operators/CracksClipper/den".
//
//    Kathleen Biagas, Fri Aug 10 13:21:06 PDT 2012
//    Support double precision data. Re-add density var to pipeline so that
//    data range can be properly set. Use pipeline var (varname).
//
//    Eric Brugger, Wed Jul 23 11:13:42 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

template <typename T1, typename T2> void
avtCracksDensityFilter_CalculateDensity(int numOriginalCells, int numNewCells, 
    T1 *emsp, T2 *volp, int* cozp, T1* denp)
{
    // find the new volume of the original cells
    double *newVol = new double[numOriginalCells];
    for (int i = 0; i < numOriginalCells; i++)
    {
        newVol[i] = 0;
    }
    for (int i = 0; i < numNewCells; i++)
    {
        newVol[cozp[i]] += volp[i];
    }

    double *volSum = new double[numNewCells];
    for (int i = 0; i < numNewCells; i++)
    { 
        volSum[i] = newVol[cozp[i]]; 
    } 

    for (int i = 0; i < numNewCells; i++)
    { 
        if (volSum[i] != 0.)
            denp[i] = emsp[i]/volSum[i];
        else 
            denp[i] = emsp[i];
    } 

    delete [] newVol;
    delete [] volSum;
}

avtDataRepresentation *
avtCracksDensityFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    if (varname.empty())
    {
        avtCallback::IssueWarning("Could not find necessary var name," 
            " in which to store density calculation.");
        return in_dr; 
    }
    vtkDataArray *vol = in_ds->GetCellData()->GetArray("ccvol");
    if (vol == NULL)
    {
        avtCallback::IssueWarning("Could not find necessary volume data," 
            " cannot calculate density. This filter must be used in"
            " conjunction with avtCracksClipper");
        return in_dr; 
    } 

    vtkIntArray *coz = vtkIntArray::SafeDownCast(
        in_ds->GetCellData()->GetArray("cracksOriginalZones"));
    if (coz == NULL)
    {
        avtCallback::IssueWarning("Could not find necessary original zones"
            " data, cannot calculate density.  This filter must be used in"
            " conjunction with avtCracksClipper");
        return in_dr; 
    }

    std::string massVar = atts.GetInMassVar();
    vtkDataArray *ems = in_ds->GetCellData()->GetArray(massVar.c_str());
    if (ems == NULL)
    {
        std::string msg = "Could not find mass varaible: " + massVar;
        msg += ".  Cannot calculate density.";
        avtCallback::IssueWarning(msg.c_str());
        return in_dr; 
    } 

    vtkIntArray *noc = vtkIntArray::SafeDownCast(
        in_ds->GetFieldData()->GetArray("originalNumCells"));
    if (noc == NULL)
    {
        avtCallback::IssueWarning("Could not find necessary number of zones"
            " data, cannot calculate density.  This filter must be used in"
            " conjunction with avtCracksClipper");
        return in_dr; 
    }

    vtkIdType numNewCells = in_ds->GetNumberOfCells();
    vtkDataArray *den = ems->NewInstance();
    den->SetName(varname.c_str());
    den->SetNumberOfComponents(1);
    den->SetNumberOfTuples(numNewCells);

    int numOriginalCells = noc->GetValue(0);

#define typeVol(massType) \
{ \
    if (vol->GetDataType() == VTK_DOUBLE) \
    { \
         avtCracksDensityFilter_CalculateDensity(numOriginalCells, numNewCells, \
             (massType*) ems->GetVoidPointer(0), \
             (double*) vol->GetVoidPointer(0), \
             (int*)coz->GetVoidPointer(0), \
             (massType*) den->GetVoidPointer(0)); \
    } \
    else \
    { \
         avtCracksDensityFilter_CalculateDensity(numOriginalCells, numNewCells, \
             (massType*) ems->GetVoidPointer(0), \
             (float*) vol->GetVoidPointer(0), \
             (int*)coz->GetVoidPointer(0), \
             (massType*) den->GetVoidPointer(0)); \
    } \
}

    if (ems->GetDataType() == VTK_DOUBLE)
    {
        typeVol(double);
    }
    else
    {
        typeVol(float);
    }

    vtkDataSet *rv = in_ds;
    rv->GetCellData()->SetScalars(den);
    rv->GetCellData()->RemoveArray("ccvol");
    rv->GetCellData()->RemoveArray(massVar.c_str());
    den->Delete();

    // Re-add the density variable to the pipeline here so that the 
    // data range can be set properly.
    double range[2] = {FLT_MAX, -FLT_MAX};
    GetDataRange(rv, range, varname.c_str(), false);
    avtDataAttributes &outAtts = GetOutput()->GetInfo().GetAttributes();
    outAtts.RemoveVariable(massVar.c_str());
    outAtts.AddVariable(varname.c_str());
    outAtts.SetActiveVariable(varname.c_str());
    outAtts.SetVariableDimension(1);
    outAtts.SetVariableType(AVT_SCALAR_VAR);
    outAtts.SetCentering(AVT_ZONECENT);
    outAtts.GetThisProcsOriginalDataExtents(varname.c_str())->Set(range);
    outAtts.GetThisProcsActualDataExtents(varname.c_str())->Set(range);
    
    avtDataRepresentation *out_dr = new avtDataRepresentation(rv,
        in_dr->GetDomain(), in_dr->GetLabel());

    return out_dr;
}
